Stability of Superfiow for Ultracold Fermions in Optical Lattices 
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Motivated by recent observations of superfluidity of ultracold fermions in optical lattices, we 
investigate the stability of superfluid flow of paired fermions in the lowest band of a strong optical 
lattice. For fillings close to one fermion per site, we show that superfiow breaks down via a dynamical 
instability leading to a transient density wave. At lower fillings, there is a distinct dynamical 
instability, where the superfluid stiffness becomes negative; this evolves, with increasing pairing 
interaction, from the fermion pair breaking instability to the well-known dynamical instability of 
lattice bosons. Our most interesting finding is the existence of a transition, over a range of fillings 
close to one fermion per site, from the fermion depairing instability to the density wave instability 
as the strength of the pairing interaction is increased. 



One of the fundamental nonequilibrium properties of a 
superfluid is the critical velocity beyond which superfiow 
breaks down. A closely related quantity is the critical 
flow momentum, Q c , defined as the maximum sustainable 
phase gradient in the superfluid. This critical momentum 
conveys information about important length scales in the 
superfluid as is easily seen for dilute quantum gases. For 
dilute bosonic superfluids, the Landau criterion tells us 
that superfluidity breaks down when the flow velocity 
exceeds the velocity of the Bogoliubov 'phonons' of the 
superfluid. The critical flow momentum is then easily 
shown to be the inverse healing length of the superfluid. 
For weakly paired fermionic superfluids the superfiow is 
limited by the small pairing gap. In this BCS regime, 
fermions depair at a critical flow momentum, which is the 
inverse Cooper pair size. As one tunes the interaction 
between fermions in a trapped Fermi gas through the 
BCS to BEC crossover, the critical momentum evolves 
from the depairing momentum of fermions to the inverse 
healing length of molecular bosons [l], 0] . 

A different route to the breakdown of superfluidity is 
through a dynamical instability as has been experimen- 
tally observed for bosonic atoms in optical lattices 
For weakly correlated bosons in a lattice, there is a crit- 
ical momentum, Q c ~ h(n/2a), where a is the optical 
lattice constant, at which superfiow breaks down as the 
effective superfluid density becomes negative 0] ■ (Hence- 
forth we will set h=l, 0=1.) For bosons at commensurate 
filling this dynamical instability has been shown [H, 0] to 
be smoothly connected to the equilibrium superfluid to 
Mott insulator transition of bosons [7( , in the sense that 
the critical momentum for the onset of this dynamical in- 
stability tends to zero as the system is tuned towards the 
Mott transition. Close to the Mott transition, Q c ~ l/£, 
where £ is the diverging correlation length for fluctuations 
near the superfluid to Mott insulator transition. Thus, 
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FIG. 1: (Color online) Critical flow momenta for phase stiff- 
ness dynamical instability, Q^, and density wave dynamical 
instability, Q p c , versus fermion filling / at various values of 
U/t in three dimensions [jj. For < 7r/2, the fermions de- 
pair at Q w Q^, but they depair at larger values of Q once 
Qt = 7r/2. Black dots indicate points where a sharp transi- 
tion from the phase stiffness to the density wave instability 



the critical flow momentum in this case conveys non- 
trivial information about an equilibrium quantum phase 
transition. 

In this Letter, we focus on the remarkable and rich 
physics associated with the breakdown of superfiow of 
paired fermions in an optical lattice. In the strong lat- 
tice potential limit, the ultracold fermion system can be 
described by a single band negative-?/ Hubbard model 
which we define on a d-dimensional hypercubic lat- 
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where cj CT creates a fermionic atom in the hyperfine spin 
state a on a lattice site i. This model is well known to 
have a superfluid ground state. The phase diagram in 
Fig. 1 summarizes the two main new results of our work, 
(i) For fillings far from one fermion per site (/ = 1), we 
show that there is a dynamical instability at a critical 
momentum where the effective superfluid density be- 
comes negative leading to the breakdown of superflow. 
For large U/t, this is related to the dynamical instability 
of lattice bosons, discussed earlier, while for small U /t it 
is closely related to depairing of fermions. (ii) We show 
that for / — 1 , there is an insulating crystal state which 
is degenerate (or nearly degenerate) with the uniform su- 
perfluid. For a range of fillings around / = 1, we show 
that this competing crystal state leads to the breakdown 
of superflow via a distinct dynamical instability at Q p c , 
which leads to a transient density wave as shown in Fig. 2. 
The critical momentum for this instability is the inverse 
density correlation length (at the relevant wavevector) in 
the uniform superfluid. 

To observe this dynamical instability in cold atom ex- 
periments, two ingredients are necessary. First, one needs 
counterpropagating laser beams, which are frequency- 
detuned from each other to generate a moving optical 
lattice. Indeed, a recent experiment using such moving 
lattices has successfully measured the critical current of 
paired fermions [l[ in a weak optical lattice. Second, one 
needs to be able to detect the density wave order, in- 
duced by the dynamical instability. This could be done 
using light scattering or noise [T3| measurements. Light 
scattering at intermediate timescales should be able to 
detect broadened Bragg spots associated with the emer- 
gent crystalline order. 

While our work is in part inspired by recent experi- 
ments demonstrating fermionic superfluidity in an opti- 
cal lattice [ll|, a broader motivation is to explore how 
competing broken symmetry states, such as charge den- 
sity waves, lead to a breakdown of superflow in uniform 
strongly correlated superfluids — this issue is of great 
importance for s yste ms ranging from high temperature 
superconductors [l2| to excitonic superfluids in quantum 
Hall bilayers 



We begin with the observation [14j that if the fermions 
only hop between nearest neighbor sites and if the density 
is fixed to be one fermion per site (/= 1), the Hubbard 
model Eq.JI]) possesses an extra pseudospin SU(2) sym- 
metry in addition to the usual SU(2) spin rotation sym- 
metry. Let us define the local version of Anderson's pseu- 
dospin operators [IH as 
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5(X) ff c L c i(T — !)• One can show that the total pseudospin 
operators J2i s i^i ± ; ^7 commute with the negative-?/ 



Hubbard Hamiltonian, where Si = +1,-1 on the two 
sublattices of the hypercubic lattice. This leads to an ex- 
act degeneracy between the superfluid and the "checker- 
board" (tt, . . . , 7r) charge density wave (CDW) states at 
/ = 1. In pseudospin language, the superfluid corresponds 
to ferromagnetic ordering of T ± while the CDW corre- 
sponds to antiferromagnetic ordering of T z . 

This degeneracy, and what lifts it, is most clearly re- 
vealed in the large U/t limit of Eq.{T]), which can be 
realized in cold atom experiments by going to a stronger 
optical lattice potential. In this limit, we can view the 
presence (absence) of a fermion pair at a site as a pseu- 
dospin up (down) configuration, and the Hamiltonian in 
Eq.([T]) maps onto a spin-a-half pseudospin model: 
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where J,j = 4tf,-/J7. At fi = 0, which corresponds to 
f = 2(Tf) + 1 = 1, and assuming Jy connects only neigh- 
boring sites with a strength J = 4t 2 /U, Eq.(2|) has SU(2) 
pseudospin symmetry. 

This pseudospin SU(2) symmetry is broken in favor 
of the superfluid for nonzero fi, which dopes the system 
away from / = 1 . It is lifted even when fi = if the pseu- 
dospin exchange Jy connects next neighbor sites, with 
a magnitude J' = At' 2 /U <C J, in which case the 'clas- 
sical' ferromagnetic state of T ± has a lower energy den- 
sity, compared to the antiferromagnetic state of T 2 , by 
an amount proportional to J'. The question we would 
like to address here is how the presence of this broken 
symmetry crystal as a competing ordered state close by 
in energy leads to an instability of the uniform flowing 
superfluid. We begin with a large U/t analysis, where 
the physics is fully described by our pseudospin model 
Eq.([2]); we then generalize our results to all values of 
U/t. For simplicity we set t' = 0, the t' ^ case is not 
substantially different. 

The evolution of the metastable current-carrying state 
of ([3|) as a function of the filling / and the flow momen- 
tum Q can be conveniently analyzed in the semiclassical 
limit of large pseudospin. Let H c e denote the classical 
Hamiltonian, obtained by replacing the pseudospin oper- 
ators in -ff xxz by angular momentum vectors Tj = TOj, 
where = (cos 9i, sin 9i cos fa, sin^ sin fa) are unit vec- 
tors and T denotes the magnitude of the angular momen- 
tum (T — 1/2 in our case). With this parametrization, 

T 2 

H c £ [O] = — Jij [cos 9i cos 9j— sin 9i sin 9j cos(fa—cj)j )} 
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For brevity, we present analytical results for / < 1; the 
physics is identical for / > 1 when J' = 0. 

Assuming the current flows in the ir-direction, we be- 
gin by extremizing the classical energy to find 9i = 9 = 
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arctan[ v //(2 - /)/(/ - 1)], ft - Q^, M = dJT(f - 
+7q), where 7 k = ± X)i=i,...,d cos(fci) and we have 
assumed that the filling /, rather than the chemical po- 
tential, is fixed. Stability of the superflow depends on 
whether this state is a local minimum of the energy func- 
tional or just a saddle point. To answer this question we 
expand H C £ to second order in fluctuations about the 
above solution and Fourier transform it to obtain: 
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where we have subtracted the energy of the uniform flow- 
ing state. Defining ?7Q(k) = (7k+Q +7k-q)/2, the three 
different stiffness coefficients in @ are given by: 



Pee{k) = 
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+ (7k + ^q(k))sin 2 (0)] 
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d 2 H ci .dJT 2 sm{29). 
H k ) = no* a j. a (7k+Q-7k-Q)-(5) 



09 



Here pee(k) is the inverse density susceptibility at 
wavevector k, p^(k) is the phase stiffness and p$$(k) 
is related to the current, transported by the mode k. 

The flowing state is a stable local minimum provided 
the following conditions are satisfied: 

pee(k), p^k) > 0, p 99 (k)p^(k) > \p 9<t) (k)\ 2 , Vk. (6) 

There are thus three distinct instabilities, which can lead 
to the decay of superflow. For reasons that will become 
clear below, the instabilities corresponding to vanishing 
pee or p^ are called dynamical, and we will refer to the 
critical flow momenta, where these occur, as Q p c and Qf 
respectively, while instability corresponding to the viola- 
tion of the last condition in Eq.© is called Landau or 
energetic instability Q and we will refer to the critical 
flow momentum for this instability as Q^ an . 

Both dynamical and Landau instabilities have the 
same physical origin, namely the superflow-carrying state 
becoming a saddle-point of the energy functional instead 
of a local minimum. In order to understand the difference 
between these instabilities, we turn to an analysis of the 
dynamics. Defining new variables pk = Tsin(#)(5#k = 
Pik + «P2k and q k = 5<f>u = gik + «<?2k, it is easy to 
show that p„k and q^ are canonically conjugate vari- 
ables which thus obey Hamiltonian quasiclassical equa- 
tions of motion. Rewritten in terms of p v and q u , the 
fluctuation Hamiltonian Eq. ((4|) has the form of a Hamil- 
tonian of a set of decoupled harmonic oscillators in two 
dimensions, labelled by k, with characteristic frequency 
wok = 2-v/ pee(k)p00(k) /Tsin(0), rotating with frequency 



Ok = 2 \p0$(k)\/Tsm(6), In this language the two dy- 
namical instabilities correspond to the point where the 
frequency of one or more of the oscillators becomes imag- 
inary, while Landau instability corresponds to the rota- 
tion frequency exceeding the oscillator frequency. Diag- 
onalizing the Hamiltonian by a canonical transformation 
we obtain the eigenfrequencies = c^ok ± Ok- Spectral 
stability requires pee{k) > 0,p<^(k) > 0. When either 
one of these become negative, the eigenmode frequencies 
acquire imaginary parts, which means that isoenergetic 
orbits in the phase space become open, i.e. the motion in 
the neighborhood of the pk = 0, q^ = saddle point be- 
comes dynamically unstable. The motion in the Landau- 
unstable regime is, however, dynamically stable in this 
linearized problem, since isoenergetic orbits in this case 
are closed. Once mode- coupling, due to nonlincarities 
or external potentials, is included the motion will most 
likely immediately become unstable. We however expect 
that at low temperatures and in two and three dimen- 
sions, where quantum corrections to the quasiclassical 
dynamics are weak [l6j |. the Landau instability in our 
system system will develop on much longer time scales 
compared to dynamical instabilities. 

After this preliminary discussion of the instabilities, 
we now calculate the dynamical phase diagram of our 
system as a function of the flow momentum Q and the 
filling /. Straightforward analysis of Eq.© gives the 
following results at T=l/2: 
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where 



A(f, d) = V(2rf - l) 2 / 2 (2 - /) 2 +8(l - /) 2 [1+(1 - f) 2 ]. 

Focussing on the dynamical instabilities, the most inter- 
esting feature of Eq. ([7]) is singular behavior of the critical 
momenta as functions of the fermion density. Defining 
/* = 1 - 1/V2d- 1, we find that for / < /*, Qt < QP, 
while for / > /*, QP < Qf. There is thus a sharp tran- 
sition, as a function of filling at /*, from an instability 
due to a vanishing phase stiffness p^ik) for < / < /* 
to an instability due to a vanishing inverse susceptibility 
pee{k) for /* < / < 1. In d = 1, we find /* = 0. The 
phase-stiffness-related instability involves all wavevectors 
of the form k = (k x ,0, . . . , 0) (for flow in the x-direction). 
The density modulation dynamical instability, by con- 
trast, occurs at a single wavevector (tt, . . . , n), since this 
instability is directly related to the checkerboard CDW 
state. For d = 3, this CDW related dynamical instability 
only exists for 1 - 1/V2 < / < 1. 

An interesting question is what is the state, towards 
which the uniform flowing superfluid becomes unsta- 
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ble, when Q£(/) < Q < Qf(f). Since the instability 
is characterized by a divergent density susceptibility at 
(tt, ...,7t), one might expect the resulting state to be 
a "flowing supersolid" [17| . in which the CDW coexists 
with superiluid flow. Indeed, one finds that such a state 
extremizes the classical energy and has lower energy than 
the uniform state for Q > Q p c . Analyzing the stability of 
such a state, however, we find that it is in turn dynami- 
cally unstable due to a negative superiluid stiffness. 

Since there does not appear to be a simple stable equi- 
librium state beyond the (%,... ,n) dynamical instabil- 
ity, we numerically integrate the semiclassical Landau- 
Lifshitz equations, ^ = {H ct ,fi\ + rjfi x {Ha,?*}, 
where { } denotes the Poisson bracket, in order to under- 
stand the fate of the system for Q > Q p . We assume a 
phenomenological damping term, r\ <C 1, to mimic weak 
dissipation arising from thermal excitations, which are 
relevant for experiment but not taken into account in our 
approach. The density order, expected to emerge beyond 
Q p c , can be characterized by the density structure factor 
S(G,t) = (4/L)E i , m ^(t)T- + ,(f)exp(-iG a ; i ). To il- 
lustrate the nature of the collective many body dynamics 
associated with the density modulation instability, we fo- 
cus on two physical observables: the density order param- 
eter, pG(t) = \J S(G, t)/L, at the wavevector G where 
S(G, t) is the largest at that time, and the fermion pair 
current I s (t) = (4/L) EiP?(W+i (*) " T? (t)Tf +1 (t)}. 
We observe (see Fig. 2) that starting from the uniform 
superiluid, there is an initial growth of a density modu- 
lation order parameter, with a growth rate proportional 
to Q — Q p and independent of the weak dissipation. The 
density modulation wavevector G fluctuates around it, 
which is the ordering wavevector of the CDW crystal. 
The growth of this CDW order leads to a drop in the 
flow momentum as the current decays via phase slips 
in regions where the local density reaches / = 1, giv- 
ing rise to a transient insulating CDW patch. At these 
intermediate timescales, the CDW order is unstable as 
discussed and exhibits fluctuations over time scales 1 / J 
(arising from excited collective modes) as seen in Fig. 1. 
In the presence of dissipation, the system eventually set- 
tles down into a new (lower energy) steady state which, 
typically, is a uniformly flowing state, where the CDW 
order has decayed away and the final flow momentum is 
Q < Q p . The main message from our numerics is that 
even when the "flowing supersolid" state is dynamically 
unstable, a transient density wave nevertheless emerges 
when Q > Q p . For modest dissipation, its fluctuations 
persist over a time scale set by the inverse dissipation 
rate as seen from the Fig. 2. 

We now generalize the large U/t results, discussed so 
far, to all values of U/t. For small fillings / < /* and 
large U/t we see that the first dynamical instability one 
encounters is at Qf — tt/2. However, for weaker cou- 
plings, one can view the Anderson pseudospin model as 
being applicable at length scales larger than the BCS co- 
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FIG. 2: (Color online) Evolution of the current, I s , and 
the scaled CDW order parameter, 2 x pc, for a one dimen- 
sional lattice of 200 sites with / = 0.8 (corresponding to 
Qc w 0.125-7T), Q = 0.13-7T and indicated values of dissipation. 



herence length, which is proportional to the fermion pair 
size £P air . Thus, this phase stiffness instability will occur 
at Qf ~ ir/2£? air instead of tt/2, where £ pair > 1 at weak 
coupling. From this result it is clear that, at weak cou- 
pling, this dynamical instability is simply related to the 
fermion depairing instability — once the flow momentum 
exceeds the inverse fermion pair size, the pairs will break 
up. Thus, when / < /*, there is a long wavelength dy- 
namical instability associated with vanishing p^, which 
smoothly evolves from the vicinity of the pair breaking in- 
stability Qf ~ < 1 when U/t<l, to Qt = ir/2 
when U/t ^> 1. At fermion densities / > /* the situation 
is completely different. In this case, for small U/t, there 
is still the above dynamical instability, which will occur 
close to the depairing instability. However, at larger U/t, 
this instability is preempted by the distinct density mod- 
ulation dynamical instability at Q p , which is set by the 
inverse correlation length of the incipient CDW order in 
the uniform superfluid. The transition between the de- 
pairing and the density wave instabilities occurs when 
QP^pmr ^ j e wnen t ne fermion pair size becomes 
smaller than the correlation length of the CDW order. 
This correlation length diverges at / = 1. Since these 
two dynamical instabilities involve completely different 
wavevectors and different susceptibilities, they can not be 
smoothly connected; thus for fillings /* < / < 1, there 
will be a sharp change in the character of the dynam- 
ical instability with increasing U/t. This sharp change 
should be contrasted with a smooth evolution of the 
ground state properties with increasing U/t (BCS-BEC 
crossover). Our result implies that while the evolution 
of the ground state is smooth, the evolution of the ex- 
citation spectrum in the presence of competing ordered 
states nearby in energy is not smooth, which is revealed 
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in the way the superfluid flow breaks down once the crit- 
ical phase gradient is exceeded. We have confirmed the 
above heuristic arguments by generalizing the BCS mean 
field theory to include flow and density wave order pa- 
rameter. The results of our calculations are summarized 



in Fig. [TJ details will be presented elsewhere [18 1. 
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